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1.  INTRODUCTION 

The  experimental  setting  and  technical  background  are  contained  in 
references  [1,  2].  It  is  presumed  that  the  reader  is  familiar  with  the 
background.  In  short,  consider  three  short  baseline  arrays  located, 
approximately,  at  the  vertices  of  an  equilateral  triangle.  Each  array  has  four 
hydrophone  receivers,  called  the  x,  y,  z,  and  c  phones,  placed  at  four  of  the 
corners  of  a  cube  and  forming  a  (local)  Cartesian  coordinate  system.  The  x,  y, 
and  z  array  arms  (each  having  length  D)  are  viewed  as  coordinate  axes 
stemming  from  an  origin  at  the  c  phone.  The  long  baseline  technique  uses 
one  phone  from  each  array.  For  convenience  of  presentation  we  will  use  only 
the  c  phones,  but  conceptually,  it  matters  not  which  phone  is  selected  from  a 
particular  array. 

A  single  sound  source  (ping)  is  received  by  each  of  the  three  c  phones  and, 
because  of  the  synchronous  nature  of  the  data  collection  system,  we  have 
three  sound  ray  transit  times;  call  them  t],  ti,  Also  the  three  position 
locations  (poslocs)  are  available  from  the  short  baseline  processing;  call  them 
x„  y,  and  z,  for  j=1,2,3,  respectively.  These  are  used  to  initialize  the  algorithm. 
All  three  of  these  locations  are  in  the  range  coordinate  system.  I.e.  the  two 
horizontal  components  are  measured  from  a  common  origin  and  the  vertical 
(z-value)  is  measured  from  the  water  surface,  positive  downward. 

Similarly  we  have  the  three  entrance  angles,  6],  62,  and  03,  each 
representing  the  estimated  elevation  angle  of  the  sound  ray  at  the  c  phones 
and  measured  from  the  horizontal.  The  overall  goal  is  to  modify  these 
entrance  angles  in  a  way  that  will  produce  a  common  vector  for  the  poslocs  of 
the  three  arrays.  (For  the  present  we  are  ignoring  all  other  sources  of  error 
and  focusing  on  the  question  of  array  tilt.)  The  difference  between  the  final 
and  initial  entrance  angle  at  a  given  c  phone  provides  an  array  tilt  correction 
for  the  azimuthal  direction  of  the  sound  source.  Two  such  corrections  for  two 
different  azimuths  can  be  used  to  correct  the  array's  two  orientation  angles, 
Xtilt  and  Vtilt-  h'  addition,  the  azimuth  error  can  be  used  to  adjust  Zrot- 


The  entrance  angle  and  the  transit  time  serve  as  the  customary  inputs  to 
ray  tracing  methods  that  return  the  horizontal  distance  and  the  vertical 
displacement  of  the  sound  source  measured  from  the  receiver  location.  (Ray 
tracing  software  of  both  the  isospeed  and  isogradient  type  are  contained  in 
reference  [2].)  These  methods  can  also  be  used  in  inverse  fashion.  E.g. 
reference  [2]  contains  a  program  RAYFIT  which  takes  the  positions  of  the 
source  and  receiver  as  input  and  returns  the  transit  time  and  the  entrance 
angle  as  output.  In  the  present  task  we  use  transit  time,  receiver  location,  and 
source  depth  as  input  and  compute  the  entrance  angle  and  horizontal  range 
of  the  source  as  output. 

The  locations  of  the  three  c  phones  are  known  in  range  coordinates  and 
from  these  values  we  can  compute  h],  h2,  and  h^,  the  horizontal  distances  of 
each  phone  from  the  position  of  the  sound  source  as  determined  originally  by 
the  short  baseline  methodology.  This  step  is  necessary  because  ray  tracing 
algorithms  are  two  dimensional  (range  and  depth). 

This  is  the  setting  for  the  description  of  the  algorithm  which  is  presented 
in  the  next  section.  The  results  of  testing  the  method  are  presented  in  section 
3.  Section  4  contains  some  discussion.  Fortran  source  code  is  contained  in  an 
Appendix. 

2.  METHOD 

The  algorithm  utilizes  three  dimensional  ray  tracing  in  an  iterative  way 
until  the  three  poslocs  (position  locations)  agree.  The  three  transit  times  are 
held  fixed  throughout.  The  first  step  is  to  adjust  the  entrance  angles  so  that 
the  three  vertical  values  are  the  same.  This  done  we  adjust  the  horizontal 
ranges  h\,  hi,  and  /13.  These  latter  values  are  the  radii  of  the  three  circles 
representing  the  locus  of  the  source  if  the  depth  were  correct.  See  Figure  1 .  In 
the  second  step  we  vary  the  common  depth  iteratively  until  the  three  circles 
have  a  common  intersection  point.  (If  there  is  no  common  intersection  then 
some  other  source  of  error  is  dominating.) 


3 


Case  1 


Case  2 


Figure  1 

More  specifically  our  technique  for  step  one  is  to  select  one  of  the  three 
values:  zi,  22,  and  23.;  determine  which  array  is  associated  with  that  value  and 
change  the  entrance  angles  at  the  other  two  c  phones  until  all  three  z's  are 
equal  to  this  selected  value.  Some  care  must  be  applied  in  doing  this.  It 
normally  follows  that  the  vertical  component  of  a  posloc  can  be  reduced  and 
the  horizontal  shortened  by  increasing  the  entrance  angle  at  the  phone; 
decreasing  that  angle  will  normally  increase  the  depth  and  the  horizontal. 
But  local  aberrations  in  the  sound  speed  depth  velocity  profile  could  cause 
some  mischief  and  in  a  given  instance  there  may  appear  departures  from  the 
monotone  pattern. 

Having  and  maintaining  a  common  vertical  for  the  three  poslocs  allows 
us  to  proceed  as  follows.  Above  we  have  two  versions  of  the  three  circles 
whose  radii  are  the  horizontal  ranges.  The  first  set  represents  the  case  in 
which  the  currently  determined  poslocs  are  too  shallow.  The  lowering  of  the 
common  level  (i.e.  increase  the  2  value)  will  increase  the  radii  and  this  can  be 
continued  until  the  circles  intersect  at  a  common  point.  The  second  set 
illustrates  the  case  in  which  the  common  vertical  level  is  too  deep.  In  this 
case  the  radii  must  be  decreased  (decrease  the  z  value)  until  the  circles 
converge.  In  the  former  case  the  entrance  angles  must  be  lowered  and  in  the 
latter  they  must  be  raised.  An  algorithmic  description  of  the  process  follows. 
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Algorithm 

Setting.  A  depth  velocity  sound  profile  is  available  for  all  sound  ray 
tracing  type  operations.  The  locations  of  the  three  receivers  (in  range 
coordinates)  are 


[«i(0 

A=]«2(0’ 

ko). 


for  i  =  1,2,3 


and  the  three  components  are  nominally  east,  north  and  vertical.  The  three 
initial  poslocs  for  the  sound  source  are 

for  I  =  1,2, 3 


and  the  associated  transit  times  and  entrance  angles  are 

ti  =  transit  time 
Of  =  entrance  angle 

Stepl.  Seek  a  common  vertical;  adjust  the  horizontal  ranges. 

i.  The  original  ranges  are 

hi=^l(xi-aj(i)f+(yi-a2(i)f  i  =  1,2,3. 

ii.  Compute  the  horizontal  distances  D]2,  D]3  of  the  first  array 
from  the  second  and  third  poslocs.  If  either  of  these  values  are 
greater  than  h],  then  let  Zo  =  max(zi,  Z2,  Z3).  I.e.  this  choice  is 
tantamount  to  a  presumption  of  case  1;  all  depths  need  to  be 
lowered.  Otherwise  let  zq  =  min(zi,  Z2, 23)- 

iii.  Find  the  two  arrays  for  which  z,-  ^  zq  and  execute  the  function 
NEWLCXT  for  each.  This  program  will  return  adjusted  values 
for  the  entrance  angle  and  the  horizontal  range.  These  values 
replace  their  original  counterparts  and  we  are  positioned  with 
inputs 

6i(p)  the  adjusted  entrance  angles 


for  i  =  1,2,3. 
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hi(p)  the  adjusted  horizontal  ranges 

and  all  z,  =  zq.  The  poslocs  thenrselves  are  adjusted  to 
conform  to  zq  and  the  appropriate  h^p)  with  the  use  of  their 
original  azimuths.  That  is,  the  values  z,  and  y,  are  updated.  In 
most  of  what  follows  the  argument  p  is  dropped  from  the 
notation. 

Caveat.  At  this  point  it  is  assumed  that  either  case  1  or  case  2  (see  figure)  is 
operative.  It  is  possible  that  two  circles  do  not  intersect  and  such  a 
condition  will  be  discovered  when  an  attempt  is  made  to  run  the 
function  CIRCSOLV.  If  that  program  returns  the  "discriminant 
negative"  message  then  one  must  back  up  and  vary  the  depth 
level  until  the  condition  is  corrected  before  one  can  proceed.  If 
correction  is  not  possible  then  the  input  structure  must  be 
declared  infeasible  and  the  computations  abandoned. 

Step  2.  Set  the  convergence  tolerance  £  >  0,  say  e  =  .01.  Choose  two  of  the 
circles.  Although  any  two  will  serve,  the  algorithm  selects  the  two 
circles  whose  (current)  adjusted  poslocs  are  the  closest.  The 
algorithm  then  views  these  arrays,  temporarily,  as  arrays  1  and  2. 
(The  affixed  h'z  and  fl'aO)  refer  to  the  horizontal  range  and  vertical 
component  of  location  of  the  third  array  in  this  temporary  index 
system.)  Solve  for  the  two  points  that  they  have  in  common.  The 
function  CIRCSOLV  will  do  this,  returning  {uj,  Vj)  for  j  =  1,  2.  The 
intersection  that  we  want  is  the  one  that  is  closest  to  the 
horizontal  location  of  the  third  array.  Thus  let 

Di=V(wi-fl'i(3))2  +  (ui-«'2(3))2 

Di  =  V(W2  'i(3))2  +  (U2  -  a'2(3))2 


D3  =  min(Di,  D2). 

Next  determine  whether  the  three  circles  have  a  common 
intersection,  or  whether  case  I  or  case  2  is  operative. 

Ut 
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SIGN=  1  ifD3>h'3  +  e 
=  -1  if  D3  <  /I'a  -  e 

=  0  if  iDa-ft'al  <£ 

Notice  that  we  have  case  1  if  SIGN  =  1,  case  2  if  SIGN  =  -1,  and  the 
three  poslocs  agree  if  SIGN  =  0.  If  SIGN  =  0  we  exit  the  program. 
In  so  doing  a  new  set  of  azimuths  will  be  determined  from  the 
common  intersection  point.  If  SIGN  0,  set  zqo  ^  zo  arid  update 
zq  with  the  equation 

20  <-  200  +  SIGN  •  I  l-h'3/D3 1  (fl'3(3)-2oo)  •  GAM 

and  GAM  is  a  tuning  constant  to  speed  convergence.  (For  our 
applications  GAM  =  5  appears  to  work  well.) 

Having  a  new  value  for  20  we  must  now  execute  NEWLOC  for  all 
three  arrays.  Thus  use  the  adjusted  horizontal  ranges  and  the 
original  azimuths  to  compute  the  adjusted  poslocs;  return  to  the 
beginning  of  Step  2. 

Supporting  functions 

NEWLOC(0,  t,  h,  2, 20,  dip),  hip)) 

Inputs  9  =  elevation  angle 

t  =  transit  time 

h  =  horizontal  range 

2  =  depth  of  posloc 

20  =  depth  goal 

Outputs 

6ip)  =  adjusted  elevation  angle 
hip)  =  adjusted  horizontal  range 

Algorithm 

Make  the  trial  adjustment 
dip)  <-  0+  izo-z)/h. 
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Call  a  ray  tracing  program  that  converts  the  trial 
entrance  angle  B(p)  and  the  fixed  transit  time  t  into 
new  values  for  h  and  z. 

Return  to  *  if  |  z  -  zq  1  >£. 

Otherwise  set  h(p)  =  h  and  exit  because  of 
convergence. 


QUAD(A,  B,  C,  Xi,  X2,  ER)  a  quadradic  equation  solver. 

Inputs  A,B/C:  the  coeficients  of  the  quadratic  equation. 

Outputs  Xi,  X2:  the  real  solutions, 

D<r-  B2-4AC 
Flag  ER  if  D  <  0 
Xi  =  (-B  +^[D)/2A 
X2=(-B  -yfD)/2A 


CIRCSOLVfA],  A2,  h],  hi,  u\.  t>i.  uj. 

faid)] 

Inputs  ^i=|ai(2)f 


h\,  hi 


V2) 

east  coordinates  of  the  two  circles 

north  coordinates  of  the  two  circles 
circle  radii 


Outputs  (w’,,  ui)  first  intersection  point 

(W2/  t>2)  second  intersection  point. 
Code  Bq  =  ir-  h\-h\  +  d\{2)  -  a^Cl)  +  02(2)  -  ^liV) 

Bi  <-2(ai{2)-ai(l)) 

B2<-2(fl2(2)-a2(l)) 

Bi  < — B1/B2  ;  Bo<-Bo/B2 
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Let 


Pi  =  1  + 

P2  =  2Bi(Bo-a2(l))-2fli(l) 
P3  =  a](l)  +  (Bo-a2il))^-h] 


Call  QUAD 


(P\,  Pi,  Ps-  u\,  U2j  ERR) 


Then  set 


end. 


U2<r-  Bo  +  BlUl 

i>2  <-  Bo  +  Bi  i>i 


Algebraic  description  of  the  solution  for  the  intersection  of  two  circles. 
Given  the  two  circles 

[x-fli(l)]2  +  [i/-<J2(1)]2  =  /J? 

U-iZl(2)]2  +  [y-^2(2)]2  =  ;72 

We  begin  by  expanding  the  left  hand  sides  and  then  subtract  the  second 
from  the  first.  This  produces 

2A:[ai(2)  - fljd)]  +  2y[a2(2)  - a2(l)]  =h]-h\^  [fli(2)  - fli(l)]  +  [al{2)  - ^2(1)] 


or 


BiX  +  Biy  =  Bq 

and  Bo,  Bi,  and  B2  are  defined  by  identification  above.  If  we  redefine 

Bi  < - B1/B2  and  Bo  <—  Bq/B2  then  the  linear  relationship  immediatelv  above 

may  be  re-expressed  as 

\/  =  Bo  +  Bjx 

and  this  may  be  substituted  for  y  in  the  original  first  circle.  Thus 
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[x-ai(l)]2  +  [Bo  +  Bix  -  =  h  j 

which  is  a  second  order  expression  in  x.  In  standard  form  it  is 

P]  +  P2X  +  P3 

where 

Pi  =  l  +  Bi 

P2  =  2Bi[Bo-a2(l)]-2fli(l) 

P3  =  4(l)+lBo-fl2(l)P-/j5 

If  the  solution  of  this  quadratic  is  real,  then  the  circles  intersect  and  the  two 
solutions  are  the  jc-components  of  those  two  vectors.  The  y-components  are 
found  from 


y  =  Bo+Bia:. 


3.  TESTING  THE  ALGORITHM 

This  section  documents  the  testing  of  the  algorithm  using  the  six  triple 
crossover  sets  available.  The  first  five  cases  are  for  10  May  1989  and  the  sixth 
for  6  June  1989.  (They  are  identified  as  cases  1,  3,  4,  5,  6,  7;  the  original  case  2  is 
defective.  The  originals  can  be  found  in  Ref.  [1],  pp.  14-17.)  Graphs  of  the  DV 
profiles  can  be  found  in  Ref.  [2]. 

The  particulars  of  testing  are  described  below.  A  filtered  central  posloc  is 
developed  (data  are  averaged)  for  each  array  in  a  triple  overlap  pointset  count. 
The  transit  time  to  that  position  is  generated  by  our  rayfitting  algorithm. 
Then  long  baseline  methods  are  used  to  produce  the  (consensus)  posloc.  We 
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also  produce  the  azimuth  and  elevation  angles  used  for  all  pjslocs.  From 
these  we  get  ZROT  corrections  and  tilt  correction  for  the  particular  azimuth. 
Azimuth  is  measured  counter  clockwise  from  the  east. 

Remark:  This  methodology  assumes  that  the  array  c-phone  locations  are 
correct.  We  used  Table  B.2  (corrected)  of  Ref.  [2].  Angles  are  measured  in 
radians  unless  otherwise  specified. 


Case  1  (5/10/89) 

ARRAY  1 

posloc 

23195.5 

-1927.8 

-326.4 

azimuth  (original) 

-.439507 

(longbase) 

-.442562 

horizontal  error  (ft) 

16.5 

elevation  (original) 

.223916 

(longbase) 

.234378 

error  (degrees) 

0.6° 

Case  3  (5/10/89) 

ARRAY  1 

posloc 

23101.0 

-1938.2 

-325.2 

azimuth  (original) 

-.451718 

(longbase) 

-.454869 

horizontal  error  (ft) 

16.8 

elevation  (original) 

.228946 

(longbase) 

.239260 

error  (degrees) 

0.6° 

ARRAY  2 

ARRAY  11 

LONGBASE 

23155.2 

23172.2 

23180.4 

-1940.2 

-1915.4 

-1934.5 

-382.1 

-364.9 

-283.4 

-2.69574 

1.56684 

-2.694418 

1.565204 

25.8 

20.8 

.209958 

.19348 

.233427 

.211358 

1.3° 

1.0° 

ARRAY  2 

ARRAY  11 

LONGBASE 

23061.2 

23076.8 

23085.6 

-1949.5 

-1928.8 

-1944.9 

-380.8 

-352.8 

-283.7 

-2.703184 

1.587809 

-2.701741 

1.585941 

24.8 

18.3 

.205557 

.196589 

.228152 

.211792 

1.3° 

o 

o 
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Case  4  (5/10/89)  ARRAY  5  ARRAY  6  ARRAY  15  LONGBASE 


posloc 

53014.7 

-2185.7 

-381.2 

azimuth  (original) 

-.552438 

(longbase) 

-.548752 

horizontal  error  (ft) 

20.3 

elevation  (original) 

.212758 

(longbase) 

.225691 

error  (degrees) 

o 

o 

Case  5  (5/10/89) 

ARRAY  16 

posloc 

64676.0 

-8633.0 

-351.0 

azimuth  (original) 

-.541822 

(longbase) 

-.548626 

horizontal  error  (ft) 

31.4 

elevation  (original) 

.202977 

(longbase) 

.210674 

error  (degrees) 

o 

o 

52981.2 

53009.1 

53011.3 

-2144.5 

-2155.8 

-2165.7 

-361.0 

-367.9 

-327.6 

-2.655442 

1.623004 

-2.649183 

1.622606 

36.8 

10.1 

.197559 

.212733 

.204975 

.222314 

0.4° 

0.5° 

ARRAY  17 

ARRAY  26 

LONGBASE 

64642.6 

64637.4 

64653.3 

-8657.7 

-8647.4 

-8654.7 

-359.5 

-343.8 

-316.8 

-2.598895 

1.529189 

-2.598187 

1.525426 

29.3 

17.5 

.216093 

.197164 

.226356 

.203495 

0.6° 

o 

o 

12 


Case  6  (5/10/89) 

ARRAY  3 

ARRAY  4 

ARRAY  54 

LONGBASE 

posloc 

38826.5 

38794.9 

38828.6 

38811.5 

915.3 

939.9 

954.0 

951.1 

-349.3 

-363.4 

-320.1 

-322.3 

azimuth  (original) 

.223151 

2.842790 

-1.393582 

(longbase) 

.231813 

2.839161 

-1.397459 

horizontal  error  (ft) 

38.8 

20.0 

16.8 

elevation  (original) 

.207200 

.269114 

.183713 

(longbase) 

.213329 

.281344 

.183226 

error  (degrees) 

0.4° 

o 

O 

-0.03° 

Case  7  (6/6/89) 

ARRAY  1 

ARRAY  2 

ARRAY  11 

LONGBASE 

posloc 

23307.2 

23265.2 

23274.3 

23288.6 

-2302.1 

-2309.4 

-2286.9 

-2304.1 

-363.5 

-392.5 

-366.2 

-298.7 

azimuth  (original) 

-.505522 

-2.607256 

1.542357 

(longbase) 

-.507964 

-2.605552 

1.538843 

horizontal  error  (ft) 

18.7 

24.0 

22.4 

elevation  (original) 

.201174 

.203625 

.211814 

(longbase) 

.216092 

.225676 

.228084 

error  (degrees) 

o 

o 

1.3° 

o 

C^ 

O 
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The  results  are  rather  dramatic,  but  cannot  be  taken  too  literally.  The  set 
of  original  poslocs  were  generated  by  the  processing  that  is  currently  in  use. 
Reference  [2]  contains  an  error  study  of  this  processing.  It  identifies  some 
sources  of  systematic  error  including  the  method  of  initializing  the  ray  tracing 
algorithm,  the  use  of  25  ft.  water  layer  increments,  and  isospeed  ray  tracing. 
The  presence  of  such  systematic  errors  can  mask  the  true  effect.  Recall  that 
the  original  poslocs  in  our  data  are  averages  of  those  poslocs  in  the  original 
triple  overlap  region;  the  transit  times  are  c.  -nputed  by  rayfitting  with  five 
feet  water  layer  increments  and  isogradient  ray  tracing. 

The  effect  is  illustrated  in  Figure  2,  which  treats  case  5  but  is 
representative  of  all  the  cases.  The  original  transit  times  for  13  point  counts 
in  the  triple  overlap  region  of  arrays  16,  17,  and  26  are  plotted  vs.  point  count. 
The  connecting  traces  serve  to  show  smoothness  in  the  original  transit  times. 
The  original  poslocs  are  rather  variable.  They  were  averaged  (as  were  the 
point  count  numbers)  to  produce  stability.  The  transit  times  fitted  to  the 
three  averages  are  marked  with  bold  gray  circles  in  the  figure.  They  appear 


Point  Count 

Figure  2.  Transit  times  to  C-Hydrophones 


14 


over  the  average  point  count  value  of  8109.3.  It  is  clear  that  they  do  not  fall 
on  the  smooth  traces.  Thus  we  have  some  systematic  error  in  our  angle  and 
other  assessments.  The  systematic  errors  identified  in  Ref.  [2]  surely 
contribute  to  this  condition. 

Remarks:  The  algorithm  seems  to  work  well.  One  should  note  that  the 
effects  of  small  errors  are  quite  remarkable.  A  one-half  degree  error  in  one  of 
the  horizontal  arrays  arms  represents  about  3  inches  (since  the  arm  is  30  feet 
long).  If  the  posloc  is  in  an  overlap  area  (say,  horizontal  range  >  4000  feet) 
then  the  corresponding  vertical  posloc  error  is  about  40  feet. 

Three  of  the  six  cases  involve  the  array  triple  (1,  2,  11).  These  cases  are 
independent,  i.e.,  the  point  count  sets  are  quite  different  for  cases  1  and  3;  and 
a  different  date  for  case  7.  It  is  encouraging  to  note  that  the  computed  tilt 
error  values  are  consistent.  This  supports  the  previously  stated  concern  about 
the  effects  of  systematic  errors. 

A  second  independent  estimate  of  the  rotation  error  obtained  from 
differing  (azimuth)  directions  may  or  may  not  be  consistent.  If  not,  it  may 
provide  information  about  the  array  locations.  But  consistency  is  necessary  to 
proceed  with  tilt  corrections. 

4.  SUMMARY 

The  algorithm  appears  to  work  quite  well.  In  all  cases  the  convergence 
was  rapid.  In  no  instance  was  the  discriminant  negative.  Indeed  a  trap  for  that 
caveat  has  not  yet  been  written  into  the  code. 

As  noted,  the  particular  results  themselves  should  not  be  taken  too 
literally.  The  systematic  errors  identified  for  the  current  processing  system 
are  the  most  likely  suspect  areas.  Uncertainty  in  the  array  locations  is  the  next 
concern.  The  use  of  filtered  poslocs  in  the  test  cases  also  may  be  contributing 
to  some  deception.  These  values  are  the  centroids  of  the  local  data  sets 
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determined  by  a  biased  processing  system.  Other  sources  of  bias  may  be 
operating,  e.g.,  drifts  in  the  depth-sound  velocity  profile. 

The  present  work  is  concerned  only  with  the  development  and  testing  of 
the  algorithm.  The  next  step  calls  for  some  sensitivity  testing.  It  is  important 
to  learn  how  the  method  is  affected  by  the  error  sources  mentioned.  There  is 
also  the  question  of  how  the  tool  should  be  used.  E.g.,  should  the  poslocs  be 
filtered  before  or  after  its  application.  Once  the  behavior  is  well  understood, 
we  can  turn  to  the  business  of  designing  an  experiment  for  calibrating  the 
array  tilts. 
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APPENDIX.  PROGRAM  LONGBASE 


PROGRAM  LONGBASE 


12/19/90 

LONGBASE.FOR  is  a  long  base  line  POSLOC  method  for  the  triple  crossover  data  from 
the  Nanoose  range. 

A1,A2,A3  are  input  vectors  providing  locations  in  3-D  for  the  three  involved  C- 
Hydrophones.  Specifically 

A  1(1)  are  the  downrange  (X)  coordinates. 

A2(I)  are  the  crossrange  (Y)  coordinates. 

_ A3(I)  are  the  depth  (Z)  coordinates.  ( Depth  is  positive  down  orientation  ). _ 


DIMENSION  A1(3),A2(3),A3(3),X(3),Y(3),Z(3),TH0(3),T(3) 
DIMENSION  H(3),AR  1  (2),AR2(2),HP(3),THP(3),PH0(3),PH  I  (3) 
DIMENSION  XX(3),YY(3),TH(X)(3)AR(4) 

COMMON  /SETl/  L(300),VEL(3{)0),V0(300),V1(300),DZ(300), 

*  LM(300) 

REAL*8  L,VEL,V0,V1.DZ,LM 

REAL*8  A 1  ,A2 A3,X,Y,Z,DZ0,AA  1  ,AA2,P1  ,P2,TH0,T,Z0,Z0(),H 
REAL*8  D1 2,D  1 3,D23,AR  1  ,AR2,D1  ,D2,D3,HP,THP,PH0,PH  1  ,TH()0 
REAL*8  VA,VB,UA,UB,Z1,Z2,Z3,TH1,EPS,XX,YY,RDH,DS,DM,GAM 

INTEGERM  DL1,DL2,DH,MARR,SIGN,AR,ERR 
CHARACTER*  15  DVNAME,OUTNAM 
CHARACTER*!  RES 

5  EPS  =  lD-3 
ERR  =  0 


C  GAM  is  a  convergence  tuning  constant.  Adjustment  may  speed  convergence. 
GAM  =  5.0D0 


C  Lx)ad  DVT  information.  Adjust  depth  to  be  positive  downward. 
WRnE(*,'(A\)')’  ENTER  VELOCITY  PROFILE: ' 
READ(*.230)DVNAME 

OPEN(UNIT=2.FILE=DVNAME,STATUS=’OLD’) 

1=  1 

10  READ(2,220,ERR=20,END=30)L(I),VEL(I) 

L(I)  =  -L(I) 

1  =  1+ \ 

GOTO  10 

20  WRITE(*,*)’  !  THERE  WAS  AN  ERROR  READING  THE  DVT  FILE!' 

STOP 

30  CLOSE(UNIT=2) 
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M  =  I-1 


C  Prepare  for  isogradient  ray  tracing. 

C  Form  the  set  of  layer  midpoints 
DO40I  =  l,M-l 

40  LM(I)  =  .5*(L(I)  +  L(I+1)) 

C  Form  depth  increments,  and  all  sound  velocity  slopes  and  intercepts. 
DO50I  =  l,M-2 
DZ(I)=LMa+l)-LM(I) 

V0(D  =  (LMa+l)*VEL(I)  -  LM(D*VEL(I+1))/DZ(1) 

V1(D  =  (VELG+l)  -  VEL(I))/DZa) 

50  CXINTINUE 

LM(M)  =  LM(M-l)  +  DZ(M-2) 


C  Enter  the  appropriate  values  for  the  data  runs. 

WRITE(*,’(A\)')'  ENTER  OUTPUT  FILE  NAME: ' 

READ(*,230)OUTNAM 

OPEN(UNIT=3,FILE=OUTNAM,STATUS='NEW') 

WRITE(*,200) 

55  WRITE(*.*)’  ’ 

C  Input  the  C-Hydrophone  location  values. 

WRITE(*,'(A\)')'  ENTER  THE  NUMBER  OF  THE  1  ST  ARRAY  :  ' 
READ(*,*)AR(1) 

CALL  ARRA  Y(AR(  1  ).A  1  (1  ),A2(1  ),A3(  1  ),ERR) 

WRITE(*,*)'  ' 

WRITE(*,’(A\)')'  ENTER  THE  NUMBER  OF  THE  2ND  ARRAY  : ' 
READ(*,*)AR(2) 

CALL  ARRA Y(AR(2),A1(2),A2(2),A3(2),ERR) 

WRITE(*,*)'  ' 

WRITE(*,’(A\)')'  ENTER  THE  NUMBER  OF  THE  3RD  ARRAY  :  ’ 
READ(*,*)AR(3) 

CALLARRAY(AR(3),A1(3),A2(3),A3(3),ERR) 

IF(ERR.EQ.1)THEN 

WRrrE(*,*)'  !  THERE  WAS  AN  ERROR  READING  THE  ARRAYS  !  ’ 
WRITE(*,*)'  !  PLEASE  TRY  AGAIN  !  ’ 

ERR  =  0 
GOTO  55 
ENDIF 
ERR  =  0 
WRITE(*,*)'  ' 

C  Input  initial  TOSLOC  values. 

57  WRITE(*,'(A\)’)'  WHICH  CASE  NUMBER  DO  YOU  WANT  TO  RUN  ?  (1-7): ' 
READ(*,*)AR(4) 

CALL  POSLOC(X,Y,Z,AR(4)) 

58  WRITE(*,270)X(1),Y(1),Z(1) 

WRITE(*,27 1  )X(2),  Y(2),Z(2) 

WRITE(*  .272)X(3),Y(3).Z(3) 

WRITEC*,*)’  ’ 

WRITE(*,’(A\)')'  ARE  THESE  VALUES  CORRECT  ?  (Y/N): ' 
READ(*,230)RES 
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WRITE(*,*)'  ’ 

IF((RES.EQ.’N’).OR.(RES.EQ.’n'))  THEN 

WRITE(*,'(A\)’)'  WOULD  YOU  LIKE  TO  ENTER  THEM  ?  (Y/N): ' 

READ(*,230)RES 

WRITE(*,*)'  ' 

IF((RES.EQ.'y').OR.(RES.EQ.’Y’))  THEN 
WRITE(*;(A\)')'  ENTER  THE  ISTPOSLOC  (X,Y,Z):  ’ 
READ(**)X(1),Y(1).Z(1) 

WRITE(*,'(A\)’)'  ENTER  THE  2ND  POSLOC  (X,Y,Z):  ’ 
READ(*,*)X(2),Y(2),Z(2) 

WRITE(*,’(A\)')'  ENTER  THE  3RD  POSLOC  (X,Y,Z): ' 
READ(*,*)X(3),Y(3),Z(3) 

GOTO  58 
ELSE 
GOTO  57 
ENDIF 
ENDIF 

C  Run  the  RA  YFIT  program  to  each  of  the  3  average  location 
C  values,  and  return  the  transit  times  T(I),  and  the  entrance  angles  THO(l). 

WRITE(*,*)’  ' 

WRITE(*,*)'  !  PROCESSING  !’ 

WRITE(*,*)'  ' 

DO  601=  1,3 


AA1=0.0D0 
AA2  =  A3(I) 

PI  =  DSQRT((X(I)-A1(I))**2  +  (Y(I)-A2(I))**2) 

P2  =  Z(I) 

IEST  =  0 

CALL  RA  YHT 1  (AA 1  AA2,P  1  ,P2,M,VEL,LM,DZ, VO, V 1 , 

*  T(I),TH0(I),TH1,IEST) 

THOOa)  =  THO(I) 

60  CONTINUE 

C  Compute  the  max  dz. 

DZO  =  MAX(DABS(Z(1)-Z(2)),DABS(Z(1)-Z(3)),DABS(Z(2)-Z(3))) 
C  Determine  horizontal  ranges  and  azimuths. 

DO  701=  1,3 

HG)  =  DSQRT((X(I)-A1(I))**2  +  (Y(I)-A2(I))**2) 

PH0(I)  =  DATAN2(Y(1)-A2(I),X(I)-A1(I)) 

70  CONTINUE 

WRITE(3,280)AR(  1  ),AR(2),AR(3)  AR(4) 

WRITE(3,250)A  1(  1  ),A2(  1 )  A3(  1 ) 

WRITE(3,25 1  )A  1  (2),A2(2)  ^3(2) 

WRITE(3,252)A  1  (3),A2(3),A3(3) 

WRITE(3,270)X(  1 ),  Y(  1  ),Z(  1 ) 
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WRITE(3,271)X(2),Y(2),Z(2) 
WRITE(3,272)X(3),Y(3),Z(3) 
WRITE(3,*)'  ’ 
WRITE(3.253)H(  1  ),H(2),H(3) 


C  Begin  step  1 . 

C  Choose  the  initial  common  depth  ZO. 

D12  =  DSQRT((X(2)-A1(1))**2  +  (Y(2)-A2(l))**2) 
D13  =  DSQRT((X(3)-A1(I))**2  +  (Y(3)-A2(l))**2) 
IF  ((D12.GE.H(l)).OR.(D13.GE.H(l)))  THEN 
ZO  =  MAX(Z(1),Z(2).Z(3)) 

ELSE 

ZO  =  MIN(Z(1);Z(2),Z(3)) 

ENDIF 

C  Determine  which  array  POSLOC  has  depth  ZO. 

DO  751=  1,3 

75  IF(Z(I).EQ.Z0)  MARR  =  I 
WRITE(3,254)DS,Z0 

GOTO  720 

80  D1  =  DSQRT((UA-A1(DH))**2  +  (VA-A2(DH))**2) 
D2  =  DSQRT((UB-A1(DH))**2  +  (VB-A2(DH))**2) 
D3  =  MIN(D1,D2) 


C  Test  for  the  confluence  of  the  three  circle  intersections. 
IF(DABS(H(DH)  -  D3).LT.EPS)  GOTO  100 


C  Perform  the  next  iteration. 

ZOO  =  ZO 

ZO  =  ZOO  +  (A3(DH)  -  Z00)*(1  -  H(DH)/D3)*GAM 
DS  =  D3  -  H(DH) 

WRITE(3,254)DS,Z0 
WRITE(3,253)H(  1  ),H(2),H(3) 

720  CONTINUE 
DO  901=  1,3 
IF  (I.EQ.MARR)  GOTO  90 

CALLNEWLOC(O.ODO,A3a),M,THO(I),T(I),H(I),Z(l),ZO) 
90  CONTINUE 

C  Step  2  control  transfer. 

IF(MARR.EQ.4)  GOTO  98 

C  Adjust  POSLOCs  for  new  horizontal  ranges. 

DO  951=  1,3 
IF(I.EQ.MARR)  THEN 
XXO)  =  X(I) 

YY(I)  =  Y(I) 

ELSE 
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XX(I)  =  Al(I)  +  H(I)*DCOS(PHO(I)) 
YY(I)  =  A2(I)  +  H(I)*DSIN(PHO(I)) 
ENDIF 

95  CONTINUE 


C  Find  the  two  closest  POSLOCs. 

D12  =  DSQRT((XX(1)  -  XX(2))**2  +  (YY(])  -  YY(2))**2) 
D13  =DSQRT((XX(1)  -  XX(3))**2  +  (YY(1)  -  YY(3))**2) 
D23  =  DSQRT((XX(2)  -  XX(3))**2  +  (YY(2)  -  YY(3))**2) 
DM  =  MIN(D12J)13,D23) 

IF(DM.EQ.D23)  THEN 
DLl  =2 
DL2  =  3 
DH=  1 

ELSEIF(DM.EQ.D1 3)  THEN 

DLl  =  1 

DL2  =  3 

DH  =  2 

ELSE 

DLl  =  1 

DL2  =  2 

DH  =  3 

ENDIF 

WRITE(3,255)DH,GAM 


C  Step  1  completed.  Set  control  for  step  2  operations  only. 
MARR  =  4 
AR1(1)  =  A1(DL1) 

AR2(1)  =  A2(DL1) 

AR1(2)  =  A1(DL2) 

AR2(2)  =  A2(DL2) 


C  Find  the  next  set  of  circles.  Transfer  to  the  test  for  confluence. 

98  CALL  CIRCSOLV(ARl,AR2,H(DLl),H(DL2),UA,VA,UB,VB) 
GOTO  80 

100  CONTINUE 

C  Finalize  after  convergence. 

IF(D3.EQ.D2)  THEN 
VA  =  VB 
UA  =  UB 
ENDIF 

WRITE(3,260)UA,  VA  ,Z0 

C  Compute  new  azimuths. 

DO  1101=  1,3 

PH0(I)  =  DATAN2(Y(I)-A2(I),X(I)-A1(1)) 

PHia)  =  DATAN2(VA-A2(I),UA-A1(I)) 

WRITE(3.26 1  )AR(I),PH0(I),PH  1  (I) 
WRITE(3,262)TH00(I),TH0(I) 
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no  CONTINUE 

CLOSE(UNIT=3) 

WRITE(*,’(A\)’)’  WOULD  YOU  LIKE  TO  RUN  ANOTHER  CASE  ?  (Y/N): 

READ(*,230)RES 
WRITE(*.*)’  ’ 

IF((RES.EQ.'Y').OR.(RES.EQ.y))  GOTO  5 
STOP 


200  FORMATC/,'  DEHNE  THE  LOCATION  OF  THE  C-HYDROPHONE  FOR  EACH 

*  'OF  THE  ARRAYS  OF/*  THE  TRIPLE  OVERLAP.  PLEASE’, 

*  '  NOTE  THAT  WE  USE  A  POSITIVE  DOWN’/’  ORIENTATION’, 

*  ’  FOR  DEPTH.’/ 

220  FORMAT(5X,D8.2,5X,D7.2) 

230  FORMAT(A) 

250  FORMAT(/10X,’A1(1)=’,F10.2.’  A2(1)=’.F10.2,’  A3(l  )=’,F10.2) 

251  FORMAT(10X,’A1(2)=’,F10.2,’  A2(2)=',F10.2,’  A3(2)=’,F10.2) 

252  FORMAT(10X,'A1(3)=’,F10.2,'  A2(3)='.F10.2,’  A3(3)=’,F10.2) 

253  FORMAT(10X,’H(1)=’,F11.4,’  H(2)=’,F1 1.4,’  H(3)=’,F1 1.4) 

254  FORMAT(/10X,’DS  =’,F11.4,’  ZO  =’,F11.4) 

255  FORMAT(/10X,'DH  =  ’,12,'  GAM  =  ’,F4.1) 

260  FORMAT(//10X,’U  =*,F10.2,’  V  =’,F10.2,’  ZO  =’, 

*  F10.2) 

261  FORMAT(/10X,'ARR  #’,I2,’  PHO  =  ’,F10.6,’  PHI  =  ’,F10.6) 

262  FORMAT(21X,'THOO  =  ',F10.6,’  THO  =  ’,F10.6) 

270  FORMAT(/10X,'X(1)  =’,F10.2,’  Y(l)  =’,F10.2,’  Z(l)  =’,F10.2) 

271  FORMATC  10X,'X(2)  =’,F10.2,’  Y(2)  =’,F10.2,’  Z(2)  =’,F10.2) 

272  FORMAT(10X,'X(3)  =’,F10.2,’  Y(3)  =’,F10.2,’  Z(3)  =’,F10.2) 

280  FORMATClOX, ’ARRAYS  ’,12,',  ’,12,’,  and  ’,12, 7X, ’CASE  ’,12) 

END 


SUBROUTINE  NEWLOC(A  1  ,A2,M,TH0,T,H0,Z,Z0) 


IINPUTS:  I 

Al: 

HORIZONTAL  COORDINATE  OF  SENSOR 

A2: 

VERTICAL  COORDINATE  OF  SENSOR,  POSITIVE  DOWN 

N: 

INDEX  OF  DEEPEST  LAYER  USED 

THO:  ELEVATION  ANGLE 

T; 

TRANSIT  TIME 

HO; 

HORIZONTAL  RANGE 

Z; 

DEPTH  OF  POSLOC 

ZO: 

DEPTH  GOAL 

OUTPUTS: 

THP: 

ADJUSTED  ELEVATION  ANGLE 

HP. 

ADJUSTED  HORIZONTAL  RANGE 
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COMMON  /SETl/  L(300),VEL(300),V0(300).V1(300),DZ(30()), 
*  LM(3(X)) 

REAL*8  L,VEL,VO,Vl,DZa.M 

REAL*8  THO,T,HO,Z,ZO,THP,HP,R,EPS,A  1  ,A2,TH  1 


C  Set  the  convergence  criteria. 

EPS  =  .OOlODO 
C  Compute  the  slant  range 
DOS  J  =  2,M 

IF((LM(J-1).LE.A2).AND.(LM(J).GT.A2))  N  -  J-1 
5  CONTINUE 

THP  =  TH0 
10  CONTINUE 

THP  =  THP  -  (Z0-Z)/H0 

CALL  ISOGRAD  1  (A  1  A2,T,THP,N,LM, VEL,V0, V 1  ,DZ,H0,Z,TH  1 ) 
IF  (DABS(Z-ZO).GE.EPS)  GOTO  10 


TH0  =  THP 


RETURN 

END 


SUBROUTINE  CIRCSQLV(ALA2,HLH2ULV1U2,V2) 


Determines  the  two  intersection  points  of  two  circles. 

INPUTS; 

A 1 :  EAST  COORDINATE  OF  THE  TWO  CIRCLES 
A2:  NORTH  COORDINATE  OF  THE  TWO  CIRCLES 
HLH2;  RADIUS 

OUTPUTS: 

U1,V1;  1ST  INTERSECTION  POINT 
U2,V2:  2ND  INTERSECTION  POINT 


DIMENSION  A1(2),A2(2) 

REAL*8  A1,A2,B0,BLH1,H2.ULU2,V1,V2,PI,P2,P3,B,B2 
INTEGER*4  ERR 


ERR  =  0 

BO  =  Hl**2  -  H2**2  +  Al(2)**2  -  Al(l)**2  +  A2(2)**2  -  A2(l)**2 
B1  =2*(A1(2)- Al(l)) 

B2  =  2*(A2(2)- A2(l)) 

BI  =-Bl/B2 
BO  =  B0/B2 
PI  = 1+B1**2 

P2  =  2*B  1  *(B0-A2(  1  ))-2*A  1(1) 

P3  =  Al(l)**2  +  (B0-A2(l))**2  -  Hl**2 
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CALL  QUAD(P1,P2,P3.U1.U2,ERR) 

IF  (ERR.EQ.1)  THEN 

WRITE(*,*)’  THERE  WAS  AN  ERROR  IN  QUAD  SOLUTION’ 

STOP 

ENDIF 

VI  =B0  +  B1*U1 
V2  =  B0  +  B1*U2 
RETURN 
END 

SUBROUTINE  QUAD(A,B,C>X1,X2,ERR) 


C  Solves  for  the  real  roots  of  the  quadratic  equation. 

REAL*8  A,B,C,D,X1,X2 
INTEGER*4  ERR 


D  =  B**2  -  4*A*C 

IF  (D.LT.O.ODO)  THEN 

ERR=  1 

GOTO  10 

ENDIF 

XI  =  (-B  +  DSQRT(D))/(2*A) 
X2  =  (-B  -  DSQRT(D))/(2*A) 
10  RETURN 
END 
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